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Abstract: The newly developed adaptive numerical dissipation control in spatially high order 
filter schemes [20] [16] for the compressible Euler and Navier-Stokes equations has been recently 
extended to the ideal and non-ideal magnetohydrodynamics (MHD) equations [17, 24, 23, 25]. 
These filter schemes are applicable to complex unsteady MHD high-speed shock/shear/turbulence 
problems. They also provide a natural and efficient way for the minimization of Div(B) nu- 
merical error [23, 24]. The adaptive numerical dissipation mechanism consists of automatic 
detection of different flow features as distinct sensors to signal the appropriate type and amount 
of numerical dissipation/filter where needed and leave the rest of the region free from numerical 
dissipation contamination. The numerical dissipation considered consists of high order linear dis- 
sipation for the suppression of high frequency oscillation and the nonlinear dissipative portion of 
high-resolution shock-capturing methods for discontinuity capturing. The applicable nonlinear 
dissipative portion of high-resolution shock-capturing methods is very general. The objective of 
this paper is to investigate the performance of three commonly used types of nonlinear numerical 
dissipation for both the ideal and non-ideal MHD. 


1. A New Method in Solving the Multi-Dimensional Conservative MHD Equations 


Consider the 3-D conservative and non-conservative (symmetrizable [6, 14]) forms of the ideal 
compressible MHD equations in Cartesian grids, 
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where the velocity vector u — (u,v,'w) T , the magnetic field vector B = (B Xy B y} B Z ) T , p is the 
density, and e is the total energy. The notation B 2 = B 2 + By + B 2 is used. The pressure is 
related to the other variables by 


P=( 7 - l)(e - \p{n- 2 + v 2 + w 2 ) - \[B 2 X + B 2 y + B 2 )). 


For plasmas and monatomic gases, 7 = 5/3. The vector on the right hand side of (1.2) is the 
non-conservative portion of the symmetrizable MHD equations and is frequently referred to in 
the literature as a source term vector. The authors prefer not to use this nomenclature since 
this is part of the symmetrizable form of the MHD and it is not a source term. 

The conservative and non-conservative forms can be written as 


U t + V - F = 0, 


U t + V • F - 5, 

respectively, where U is the state vector, F is the conservative inviscid flux vector tensor and S 
is the non-conservative portion of the equations in (1.2). The conservative and non-conservative 
forms of the non-ideal compressible MHD [5] takes the form 


U t + V * F = F*, 


U t + V - F = F v + 5, 

1 1 T 

= [ 0 divr div(u T r) + divh div((V x B) x B) — (AB — VdivB) ] . 

a (J 

The vector F v includes viscosity, resistivity, and conductivity with r being the viscous stress 
tensor, a the resistivity coefficient, and h the heat flux. 

2. Solving the Conservative System Using the Non- Conservative Eigenvectors 

Our adaptive numerical dissipation mechanism consists of automatic detection of different flow 
features as distinct sensors to signal the appropriate type and amount of numerical dissipa- 
tion/filter where needed and leave the rest of the region free from numerical dissipation con- 
tamination. An important ingredient in our method is the use of the dissipative portion of 
high-resolution shock-capturing schemes as part of the nonlinear filters for discontinuity cap- 
turing. If the dissipative portion of higher order Lax-Friedrichs or Nessyahu-Tadmor [9] type 
of shock-capturing schemes are not employed (see [19, 20] for a discussion), these nonlinear fil- 
ters usually involve the use of approximate Riemann solvers. Due to the fact that the inviscid 
fluxes of the 2-D and 3-D conservative MHD systems do not have a complete set of eigenvectors, 
we will, therefore, first discuss a new form of high-resolution shock-capturing schemes for the 
conservative MHD equations using the non-conservative eigensystem [17, 24]. 
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For non-ideal MHD, we apply the inviscid MHD base scheme twice for the viscous flux 
derivatives. There is no viscous filtering involved. In addition, without lost of generality, we will 
describe our numerical methods for the inviscid x-flux of the ideal MHD (1.1) on a uniform 
grid. The schemes to be discussed, in most part, only spell out the ^-component terms with the 
y- and ^-components omitted. Let A(U) denote the Jacobian dF/dU with the understanding 
that the present F and S are the inviscid x-component of the 3-D description above. For later 
discussion we write the non-conservative term S, in the x-direction as N(U)U X . 

Seven of the eigenvalues and eigenvectors are identical for the “conservative” Jacobian matrix 
A and the “non-conservative” Jacobian matrix ( A — N) [4, 1]. For ease of reference, we refer 
to the distinct eigenvalue (eigenvector) between the conservative and non-conservative MHD 
as the eighth eigenvalue (eigenvector). The eighth eigenvector of A of the conservative system 
associated with the degenerate zero eigenvalue can sometimes coincide with one of the other 
eigenvectors, thereby, maldng it difficult to obtain the Roe-type approximate Riemann solver 
for the multi-dimensional conservative MHD. On the other hand, the eigenvectors of the non- 
conservative Jacobian A* = (A — N ) always form a complete basis, and can be obtained from 
analytical formulas [6, 14] for 1-D or higher. A Roe-type average state was developed in [1] for 
the 8-wave 1-D conservative MHD and extended to the 3-D non-conservative MHD by Gallice 

[4]- 


In [17, 24, 23], we proposed to use eigenvectors of the non-conservative form but with the 
degenerate eigenvalue replaced by an entropy correction [7, 19] (a small parameter e that is 
scaled by the largest eigenvalue of A(U)) for the conservative form. For more than one-space 
dimension, a multi-dimension entropy correction as proposed in [19] is used for each of the 
degenerate eigenvalues in each spatial direction. Our rationale for doing this is that the incorrect 
eigenvector for the conservative form will be multiplied by an eigenvalue which is close to zero. 
Thus the effect of a “false” eigenvector will be small. Note that in the present context, the 
use of the an entropy correction is different from the standard entropy correction associated 
with expansion shocks in the Roe- type approximate solver in gas dynamics [7, 19], since the 
conservative inviscid gas dynamics equations are strictly hyperbolic. 


3. Description of High Order Filter Methods 

Our filter idea is very general and can be used in conjunction with spectral, compact and non- 
compact spatially central base schemes. Basically, the filter method consists of two steps, a 
divergence-free preserving (base scheme) step (not involving the use of approximate Riemann 
solvers or flux limiters) and a filter step (usually involving the use of approximate Riemann 
solvers and flux limiters). In order to have good shock-capturing capability and improved non- 
linear stability related to spurious high frequency oscillations, the blending of a high order non- 
linear filter and a high order linear filter was proposed in [22, 24]. The nonlinear filter consists 
of the product of an artificial compression method (ACM) indicator or wavelet sensor and the 
nonlinear dissipative portion of a high-resolution shock-capturing scheme. The high order linear 
filter consists of the product of another sensor, a tuning parameter and a high order centered 
linear dissipative operator that is compatible with the order of the base scheme being used. Here 
the extension with a modification of the gas dynamic filter approach to the MHD equations that 
minimizes the V- B numerical error is summarized. Due to space limitation, only the nonlinear 
filter using wavelets as sensors is presented here. The numerical examples considered exhibit 
highly accurate solutions without the use of high order linear filters. 

3.1. DIVERGENCE-FREE PRESERVING BASE SCHEME STEP 

The first step of the numerical method consists of a time step via a high order non-dissipative 
spatial and high order temporal base scheme operator L *. After the completion of a full time 
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step of the base scheme step, the solution is denoted by U* 

U* = L*(U n ), (3.3) 

where U n is the numerical solution vector at time level n. For example, a spatially divergence- 
free preserving eighth-order linear dissipation with the sixth-order centered base scheme to 
approximate F(U) X is written as (with the grid indices k and l for the y- and 2-directions 
suppressed) 

~ « D 06 Fj + d(Ax) 7 (D + D.) 4 U jy (3.4) 

where Dqq is the standard sixth-order accurate centered difference operator, and D+D_ is the 
standard second-order accurate centered approximation of the second derivative. The small 
parameter d is a scaled value in the range of 0.00001 to 0.01, depending on the flow problem, 
and has the sign which gives dissipation in the forward time direction. The Doe operator is 
modified at boundaries in a stable way by the so called summation- by-part (SBP) operators 
[11, 10, 22]. The linear numerical dissipation operator D + D_ is modified at the boundaries to 
be semi-bounded [15]. The base scheme step with the fourth-order classical Runge-Kutta time 
discretization takes the form 
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where the L operator is a semi-discrete form of the conservative system (1.1) or non-conservative 
system (1.2). For example, using formula (3.4) to discretize (1.1) with the flux tensor F = 
(F ) G : H ) takes the form 


=L(U)j,k,i 

= -mu F jXl + (d^Axyimmu^Ki (3 6) 

+(D k UG jA i + (dy)(Ay) 7 [(D k ) + (D k )-} 4 Uj,k,i 
+(DiUH m + (4)(A2) 7 [(A) + (A)-] 4 ^',m}> 

where Dj , D k and Di denote finite difference operators acting in the j-, k-, /-directions, respec- 
tively. Here, e.g., (Dj) oe denotes the sixth-order centered difference operator in the ^-direction. 
Similarly, (Dj)+(Dj)- denotes the second-order centered approximation of the second derivative 
in the j -direct ion. The small parameters d x , d y and d z are the same as d in (3.4). This highly 
accurate spatial base scheme is employed to numerically preserve the divergence-free condition 
of the magnetic field (to the level of round-off error) for uniform Cartesian grids with periodic 
boundary conditions. 


3.2. ADAPTIVE NUMERICAL DISSIPATION FILTER STEP 

After the completion of a full time step of the divergence-free preserving base scheme step, the 
second step is to adaptively filter the solution by the product of a “wavelet sensor” and the 
“nonlinear dissipative portion of a high-resolution shock-capturing scheme”. If necessary, 
the blending of a high order linear filter with a nonlinear filter [22] can be employed and not 
discuss here. The final update step after the filter step can be written as (with some of grid 
indices suppressed for ease of illustration) 



ADAPTIVE NUMERICAL DISSIPATION CONTROL IN MHD 


5 


US = Um - f^Kli/2 - H l- 1/2] - ^f<./2 - - f^i'1/2 - W ]• (3.7) 

Here, Hj ±1 / 2 , H an ^ -^±1/2 axe the filter numerical fluxes in the x, y and z-directions, 

respectively. The x-filter numerical flux vector Hj+u 2 (with the y- and ^-components grid indices 
suppressed) is 


#/+i/2 “ #j+l/2#j+l/2> 

where i? J+1 / 2 is the matrix of right eigenvectors of the Jacobian of the non-conservative MHD 
flux vector (A,. +1 / 2 — .Vj+1/2) evaluated at the Gallice average state [ 4 ] in terms of the U* 
solution from the base scheme step ( 3 . 3 ). The notation Rj+1/2 stands for Rj+i/2,k,l and the 
subscript in Rj+1/2 indicates the average state evaluated in the x-direction of the eigenvectors 
in terms of U*. See [ 4 ] or Appendix A of [ 24 ] for the average state formula for the 3 -D non- 
conservative system ( 1 . 2 ). The H. 7+1/2 (involving the use of wavelet sensors and flux limiters) are 
also evaluated from the same average state. The dimension-by-dimension procedure of applying 
the approximate Riemann solver is adopted. 

Denote the elements of the vector H j+ 1/ 2 by h l j+ 1 / 2 ,l = 1,2, ..., 8 . The nonlinear portion of 
the filter h l j +1 / 2 ,l = 1, 2, ..., 8, has the form 

hj+1/2 — 2 ( sN ) l j +1/2(4*^ +1/2)- ( 3 - 8 ) 

Here (■S'^)^ +1 / 2 the sensor to activate the higher order nonlinear numerical dissipation filter. 
For example, (s Ar )^ +1 / 2 * s designed to be zero or near zero in regions of smooth flow and near one 
in regions with discontinuities. (s^)j +1 / 2 var i es from one grid point to another and is obtained 
from a wavelet analysis of the flow solution [ 16 ]. The function (fij+1/2 is the dissipative portion 
of the nonlinear filter for the local Zth-characteristic wave [22] in the x-direction. Note that the 
wavelet sensor can be obtained from the characteristic variables for each wave or a single sensor 
for all eight waves, based on pressure and density. Both methods were implemented but for the 
numerical test in this paper, the simpler non-characteristic sensor was employed. 

The dissipative portion of the nonlinear filter ^ +1 / 2 = 9 j+i /2 ~ ^j+ 1/2 * s the dissipative 
portion of a high order high-resolution shock-capturing scheme for the local Zth-characteristic 
wave. Here yj +1/ , 2 and ^- +1 / 2 are numerical fluxes of the uniformly high order high-resolution 
scheme and a high order central scheme for the Ith characteristic, respectively. It is noted that 
b l j+ 1/2 might not be unique since there is more than one way of obtaining 4> l j+\/2- 

Three forms of nonlinear dissipation <^- +1 / 2 are considered, namely: 

— Dissipative portion of the fifth-order WENO scheme (WEN 05 ) [8]. It can be obtained e.g., 
in the x-direction by taking the full WEN 05 scheme in the x-direction and subtracting 
D 06 Fj. 

— Dissipative portion of the a second-order MUSCL scheme [20]. 

— Dissipative portion of the Harten-Yee TVD scheme [ 20 , 24 ]. 

For example, the forms of Harten and Yee and symmetric TVD schemes [ 19 ] are already in 
the proper form in the sense that they are written in a central differencing portion ^ +1 / 2 and a 

nonlinear dissipation portion 4> l j+i/2 ■ No work is required to obtain ^ +1 / 2 in this case. 

The nonlinear filter given by ( 3 . 8 ), if applied to the entire MHD system, will not preserve the 
divergence free magnetic field condition in general. For the computations in this paper, the “No 
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filter on B” option is chosen. The nonlinear filter step of (3.8) only acts on the gas dynamic 
portion of the equations. That is, the nonlinear filter step (3.8) only applies to the first five 
equations of (1.1) or (1.2). Here the complete set of eigenvalues and eigenvectors of the full non- 
conservative MHD system is used to evaluate the first five equations of (1.1) or (1.2). With the 
divergence free spatial base scheme, the divergence free property should be preserved. Extensive 
grid convergence comparison of the “no filter on B” with the “filter all of the MHD equations” 
(filter all) options were presented in [24]. Alternative approaches in obtaining divergence- free 
preserving shock-capturing filters follow in a similar vein as the constrained transport approach 
[3] or the projection method [26]. 


4. 2-D Compressible Ideal and Non-ideal MHD Numerical Examples 

The wavelet filter schemes using the dissipative portion of WENO 5', second-order MUSCL and 
Harten-Yee TVD schemes with sixth-order spatial central base scheme (d = 0 in (3.4)) for both 
the inviscid and viscous MHD flux derivatives and a fourth-order Runge-Kutta method are de- 
noted by WAV66weno5, WAV66mus and WAV66hy respectively. The first number indicates the 
order of the base scheme for discretizing the inviscid flux derivatives. The second number indi- 
cates the order of the scheme for discretizing the viscous flux derivatives, if present. As mentioned 
before, there is no filtering for the viscous fluxes. If an eighth-order linear dissipation (d ^ 0 in 
(3.4)) is used for the base scheme, the symbol “AD8” is added as in WAV66weno5+AD8. Com- 
putation using the same temporal and spatial scheme for the viscous MHD flux derivatives, and 
the standard fifth-order WENO scheme for the inviscid flux derivatives is denoted by WEN05. 
Computations using a second-order MUSCL and the Harten-Yee [24] TVD scheme for the in- 
viscid MHD flux with the second-order central scheme for the viscous flux and a second-order 
Runge-Kutta method axe denoted by MUSCL and HY, respectively. 

The entropy fix parameter e is 0.25 for the Harten-Yee, MUSCL, WAV66mus and WAV66hy 
schemes (to avoid expansion shocks and carbuncle phenomenon) . The cut off wavelet Lipschitz 
exponent j3 is 0.5 [16] for all the wavelet filter schemes. See [20, 16, 22] or Appendix B of [24] 
for the definition of e and /?. Except for WEN05, the van Leer version of the van Albada limiter 
is used. For the second-order MUSCL scheme, the limiter is applied to the primitive variables. 
For viscous computations, the Reynolds number Re = 1000, resistivity coefficient of 100 and 
a Prantdl number of 0.72 are used for all test cases. Extensive grid convergence studies using 
WAV66hy and ACM66hy (using ACM instead of wavelet as the sensor) for typical ideal and 
non-ideal MHD test cases were conducted in [17, 24, 25]. More accurate solutions were obtained 
with WAV66hy and ACM66hy than with WEN05, which is more CPU intensive. The following 
investigates the performance of the three different filters. 

4.1. COMPRESSIBLE MHD ORSZAG-TANG VORTEX (7 = 5/3, PERIODIC BC) 

The 2-D Compressible MHD Orszag-Tang vortex problem [12, 13, 2] consists of periodic bound- 
ary conditions with smooth initial data . 

(p, u, v, w, p, B x , By, B z ) = (25/9, — siny, sinx, 0, 5/3, — siny, sin2x, 0). 

The initial sine waves break into discontinuities at a later time with complicated flow interactions. 
The computational domain is 0 < x < 2n, 0 < y < 27r and the computation stops at time T= 
3.14 (~ 7r), when discontinuities have formed and interacted. The solution has both complicated 
structure and discontinuities. Density contours with 30 equally spaced contours between 0.9 and 
6.1 are used for illustration. 

Figures 1 and 2 show the comparison among the three filter schemes (no filter on B option), 
WENO 5, MUSCL and Harten-Yee (HY) using uniform 101 x 101 and 801 x 801 grids for the 
inviscid MHD. Figures 3 and 4 show the corresponding viscous MHD comparison. 
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MUSCL, 101x101 Harten-Yee, 101x101 



* x x 
Fig. 1. Inviscid Orszag-Tang Problem using a 101 x 101 grid. Density contours for the 

ideal MHD at T = 3.14. Top row: WEN05(left), MUSCL (middle), HY (right). Bottom row: 

WAV66weno5(left), WAV66mus (middle), WAV66hy (right). 



X x * 

Fig. 2. Inviscid Orszag-Tang Problem using a 801 x 801 grid . Density contours for the 
ideal MHD at T = 3.14, Top row: WENOSfleft), MUSCL (middle), HY(right). Bottom row: 
WAV 66weno5 (left), WAV 6 6mus (middle), WAV 6 6hy (right). 


Convergence solutions were obtained by all six methods using’ the 801 x 801 grid. Computa- 
tions based on a 1601 x 1601 grids are used as the reference solutions. For 51 x 51 and 101 x 101 
grids, small structures are more well captured by the three filter methods than that of WEN05, 
MUSCL and Harten-Yee. In addition, for the inviscid case, the three filter methods are more 
stable than that of the other three methods in the sense larger CFL number can be used. Figure 
2 shows the computations using a CFL of 0.6 and a 801 x 801 grid. WEN05 and MUSCL show 
a slight small oscillation. By applying the limiter in the characteristic variables in the MUSCL 
scheme these oscillations are suppressed. 
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MUSCL, 101x 101, viscous Harter>-Yee, 101x 101, viscous 



x x x 
Fig. 3. Viscous Orszag-Tang Problem using a 101 x 101 grid. Density contours for 

non-ideal MHD with Re = 1000, and resistivity coeff. of 100 at T = 3.14. 

Top row: WEN05(left), MUSCL (middle), HY(right). Bottom row: WAV66weno5(left), 

WAV 6 6mus (middle), WAV 6 6hy (right). 
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X x x 

Fig. f. Viscous Orszag-Tang Problem using a 801 x 801 grid. Density contours for 
non-ideal MHD with Re — 1000, and resistivity coeff. of 100 at T = 3.14. 

Top row: WEN05(left), MUSCL (middle), HY (right). Bottom row : WAV66weno5(left), 
WAV66mus (middle), WAV66hy (right). 

For the viscous case, the flow structure is less complicated than that of the inviscid case. All 
computations use a CFL of 0.6. For coarse grids, again small structures are better captured by 
the three filter methods than by WEN05, MUSCL and Harten-Yee. For both the inviscid and 
viscous computations, all three filter methods using the no filter on B option are divergence- 
free preserving. Whereas the “filter all” option as well as WEN05, MUSCL and HY are not 
divergence free. Their V- B numerical error at T = 3.14 increases as the grid is refined. 
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For this test case, MUSCL and Harten-Yee require similar CPU time. The CPU required 
by the three filter methods are within 15% depending on the problem, grid spacings and time 
steps. They require slightly more CPU time (20%) than the Harten-Yee and MUSCL schemes. 
WEN05 requires at least twice the CPU time of all considered methods. 

4.2. A PLANAR SHOCK INTERACTING WITH A MAGNETIC CLOUD (7 = 5/3, SUPERSONIC 
INFLOW & OPEN BOUNDARIES) 

The second test problem is a planar shock interacting with a magnetic cloud studied in [ 2 ] . This 
is a more challenging problem to simulate due to the stiffness of the flow with rapidly developing 
complex wave interactions close to the boundaries. The same initial configuration as in [18] is 
considered here. The computational domain is the square 0<x<l,0<y<l. A planar shock 
is initially situated at x = 0.6 and moves towards the right. 


(p, U, V , w, p, B x , By, B z ) l = (3.86859, 0, 0, 0, 167.345, 0, 2.1826182, -2.1826182), 

(p, u, V , to, p, B x , By, B z ) r = ( 1 , -11.2536, 0, 0, 1, 0, 0.56418958, -0.56418958). 

Inside the cloud with center at (0.8, 0.5) and radius 0.15, a state of increased density is given 
by the initial state: 


(p, u, v, to, p, B x , B y , B z ) c = (10; -11.2536, 0, 0, 1, 0, 0.56418958, -0.56418958). 

The flow velocity is directed in the negative x-direction, and the cloud will move to the left. The 
right boundary is supersonic inflow, where the right state is imposed. The other boundaries are 
open boundaries. Density contours with 50 equidistant contours in log scale from log(0.99) to 
log(48) are used. 

Figures 5 and 6 show the comparison among the three filter schemes (WAV66weno5+AD8, 
WAV 66 mus+AD 8 , WAV 66 hy-|-AD 8 , d=0.001 in (3.4)), with the “no filter on B option”, WEN05, 
MUSCL and Harten-Yee using uniform 201 x 201 and 801 x 801 grids for the inviscid MHD. 
Computations using WAV66weno5, WAV 66 mus, and WAV 66 hy (d=0 in (3.4)) are not as sta- 
ble. Figures 7 and 8 show the corresponding viscous MHD comparison using WAV66weno5, 
WAV 66 mus, and WAV 66 hy (d=0 in (3.4)). For the viscous case, AD 8 (i.e. , d 0) is not neces- 
sary as part of the base scheme for a stable solution. 

As opposed to the Orszag-Tang problem, in this case, grid convergence was not quite achieved 
for all six methods using the 801 x 801 grid for the inviscid. For both the viscous and inviscid 
cases, overall the three filter methods are more accurate than the other three methods. Compu- 
tations based on a 1601 x 1601 grid are used as the reference solutions (see Fig. 9). For 51 x 51 and 
101 x 101 grids, small structures are better captured by the three filter methods than WEN05, 
MUSCL and Harten-Yee. Moreover, among the three filter methods, WAV 66 hy+AD 8 is slightly 
more accurate on capturing the fine scale structures in most coarse grid cases. It is interesting 
to see the contrast in the flow structure between the inviscid and viscous converged solutions. 

For all three filter methods using the no filter on B option, perfect V • B preservation within 
machine zero is only obtained up to a certain time (T = 0.04). The increase in the norm of V - B 
is caused by boundary effects. A SBP difference boundary operator is used [11, 10]. Due to the 
wide grid stencil of the SBP boundary difference operator in conjunction with the need to use 
an extrapolation to the outermost open boundary point, V- B is not preserved. The effect is 
only seen when the solution is non-trivia! on the boundary where complex wave interactions axe 
taking place in both directions of the open boundaries. Although divergence- free preservation by 
the three filter methods is not possible for T > 0.04 by the “no filter on B” option, the L 2 -norm 
of V - B for this option is at least an order of magnitude smaller than the “filter all” option and 
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X X X 

Fig. 5. Inviscid Magnetic Cloud Problem using a 201 x 201 grid. Density contours for the 
ideal MHD at T — 0.06. Top row: WEN05(left), M US CL (middle), HY(right). Bottom row: 
WAV66weno5+AD8(left), WAV66mus+AD8(middle), WA V66hy+AD8( right) . 



X X X 

Fig. 6. Inviscid Magnetic Cloud Problem using a 801 x 801 grid. Density contours for the 
ideal MHD at T = 0.06. Top row : WENOS(left), M US CL (middle), HY(right). Bottom row: 
WAV66weno5pAD8(left), WAV66mus+AD8 (middle), WAV66hy+AD8(right). 

the three standard shock-capturing schemes when T > 0.04. See [24, 23] for some illustrations 
of WAV66hy+ AD8. An alternative in obtaining a divergence-free solution is to employ standard 
divergence free cleaning on the filter step. This is a subject of on going research. The relative 
CPU required by the six methods are the same as the previous test case. 

5. Concluding Remarks 

The performance of three nonlinear filters for the adaptive numerical dissipation control in high 
order methods [20, 16, 17, 24] is investigated. The “no filter on B” option by the three filter 
methods works well for both the conservative and non-conservative (computations not shown) 
systems and exhibits smaller V* B numerical error than standard shock-capturing methods with- 
out traditional divergence cleanings. For periodic boundary conditions and for open boundaries 
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Fig. 7. Viscous Magnetic Cloud Problem using a 201 x 201 grid. Density contours 
for non-ideal MHD with Re = 1000, and resistivity coeff. of 100 at T = 0.06. 

Top row: WEN05(left), MU S CL (middle), HY(right). Bottom row: WAV66weno5(left), 

WAV66mus (middle ), WAV66hy (right). 



XXX 

Fig. 8. Viscous Magnetic Cloud Problem using a 801 x 801 grid. Density contours 
for non-ideal MHD with Re = 1000, and resistivity coeff. of 100 at T — 0.06. 

Top row: WEN05(left), MUS CL (middle), HY(right). Bottom row: WAV66weno5(left), 

WAV 6 6mus (middle), WAV66hy (right). 
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Fig. 9 . Inviscid and Viscous Magnetic Cloud Problem using a 1601 x 1601 grid by 
WAV66weno5+AD8 ( d=0.001 ). Density contours for ideal and non-ideal MHD with Re = 1000, 
and resistivity coeff. of 100 at T — 0.06 ; inviscid (left) and viscous (right) 
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without complex wave interactions near the physical boundaries, these filter schemes are diver- 
gence free. In general, for coarse grids, the high order methods are more accurate and require 
less grid points than required by second-order methods. For fine enough grids, in most test cases, 
the accuracy is similar for all six methods. Among the three filter methods, WAV66hy+AD8 is 
slightly more accurate on capturing the fine scale structures in most of the coarse grid cases. 

For the two test cases, MUSCL and Harten-Yee require similar CPU time. The CPU required 
by the three filter methods are within 15% depending on the problem, grid spacings and time 
steps. They require slightly more CPU time (20%) than the Harten-Yee and MUSCL schemes. 
WEN05 requires at least twice the CPU time of all other of the methods. This is due to the fact 
that all filter schemes require only one Riemann solve per time step per direction (independent of 
the time discretizations of the base scheme step) as opposed to two Riemann solves per time step 
per direction by the MUSCL, Harten-Yee schemes using a second-order Runge-Kutta method. 
In addition, for the two test cases and all six methods (except the no filter on B option for the 
three filter schemes), the V* B numerical errors (at their corresponding stopping times) increase 
as the grid is refined. 
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